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Abstract 

Secular variations of the geomagnetic field have been measured with a continuously im- 
proving accuracy during the last few hundred years, culminating nowadays with satellite data. 
It is however well known that the dynamics of the magnetic field is linked to that of the velocity 
field in the core and any attempt to model secular variations will involve a coupled dynamical 
system for magnetic field and core velocity. Unfortunately, there is no direct observation 
of the velocity. Independently of the exact nature of the above-mentioned coupled system 
- some version being currently under construction - the question is debated in this paper 
whether good knowledge of the magnetic field can be translated into good knowledge of core 
dynamics. Furthermore, what will be the impact of the most recent and precise geomagnetic 
data on our knowledge of the geomagnetic field of the past and future? These questions are 
cast into the language of variational data assimilation, while the dynamical system considered 
in this paper consists in a set of two oversimplified one-dimensional equations for magnetic 
and velocity fields. This toy model retains important features inherited from the induction 
and Navier-Stokes equations: non-linear magnetic and momentum terms are present and its 
linear response to small disturbances contains Alfven waves. It is concluded that variational 
data assimilation is indeed appropriate in principle, even though the velocity field remains 
hidden at all times; it allows us to recover the entire evolution of both fields from partial and 
irregularly distributed information on the magnetic field. This work constitutes a first step on 
the way toward the reassimilation of historical geomagnetic data and geomagnetic forecast. 



1 Introduction 



The magnetic observation of the earth with sateUites has now matured to a point where continuous 
measurements of the field are available from 1999 onwards, thanks to the Oersted, SAC-C, and 
CHAMP missions (e.g. Olsen et al. 2000 Maus et al. 2005 and references therein) . In conjunction 



with ground-based measurements, such data have been used to produce a main field model of 



remarkable accuracy, in particular concerning the geomagnetic secular variation (GSV)( Olsen 
et al. 2006a I . Let us stress that we are concerned in this paper with recent changes in the earth's 



magnetic field, occurring over time scales on the order of decades to centuries. This time scale is 
nothing compared to the age of the earth's dynamo (> 3 Gyr), or the average period at which the 



dynamo reverses its polarity (a few hundreds of kyr, see for instance iMerrill et al. (1996)), or even 



the magnetic diffusion time scale in earth's core, on the order of 10 kyr (e.g. Backus et al. 19961 



It is, however, over this minuscule time window that the magnetic field and its changes are by far 

T989J . 



best documented (e.g. ^Bloxham et al 
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Downward-projecting the surface magnetic field at the core-mantle boundary, and applying 
the continuity of the normal component of the field across this boundary, one obtains a map of 
this particular component at the top of the core. The catalog of these maps at different epochs 
constitutes most of the data we have at hand to estimate the core state. Until now, this data has 
been exploited within a kinematic framework (Roberts and Scott 1965 Backus 1968 ): the normal 
component of the magnetic field is a passive tracer, the variations of which are used to infer the 

For the purpose of modeling the 



velocity that transports it (e.g. Le Mouel 



Le Mouel 


1984 


Bloxham 


1989 



core field and interpreting its temporal variations not only in terms of core kinematics, but more 
importantly in terms of core dynamics, it is crucial to make the best use of the new wealth of 
satellite data that will become available to the geomagnetic community, especially with the launch 



of the SWARM mission around 2010 (Olsen et al. 2006b I 



This best use can be achieved in the framework of data assimilation. In this respect, geomag- 
netists are facing challenges similar to the ones oceanographers were dealing with in the early 
Nineteen-nineties, with the advent of operational satellite observation of the oceans. Inasmuch 
as oceanographers benefited from the pioneering work of their atmosphericist colleagues (data 
assimilation is routinely used to improve weather forecasts), geomagnetists must rely on the de- 
velopments achieved by the oceanic and atmospheric communities to assemble the first bricks of 
geomagnetic data assimilation. Dynamically speaking, the earth's core is closer to the oceans than 
to the atmosphere. The similarity is limited though, since the core is a conducting fluid whose 
dynamics are affected by the interaction of the velocity field with the magnetic field it sustains. 
These considerations, and their implications concerning the applicability of sophisticated ocean 
data assimilation strategies to the earth's core, will have to be addressed in the future. Today, 
geomagnetic data assimilation is still in its infancy (see below for a review of the efforts pursued 
in the past couple of years). We thus have to ask ourselves zero-th order questions, such as: 
variational or sequential assimilation? 

In short, one might be naively tempted to say that variational data assimilation (VDA) is 
more versatile than sequential data assimilation (SDA), at the expense of a more involved im- 
plementation -for an enlightening introduction to the topic, see Talagrand (19971. Through an 
appropriately defined misfit function, VDA can in principle answer any question of interest, pro- 
vided that one resorts to the appropriate adjoint model. In this paper, we specifically address 
the issue of improving initial conditions to better explain a data record, and show how this can 
be achieved, working with a non-linear, one-dimensional magneto-hydrodynamic (MHD) model. 
SDA is more practical, specifically geared towards better forecasts of the model state, for example 
in numerical weather prediction (Talagrand 1997). No adjoint model is needed here; the main 



difficulty lies in the computational burden of propagating the error covariance matrix needed to 
perform the so-called analysis, the operation by which past information is taken into account in 



order to better forecast future model states (e.g. Brasseur 2006 ) 



Promising efforts in applying SDA concepts and techniques to geomagnetism have recently been 
pursued: |Liu et aL (20071 have performed so-called Observing System Simulation Experiments 
(OSSEs) using a three-dimensional model of the geodynamo, to study in particular the response 
(as a function of depth) of the core to surface measurements of the normal component of the 
magnetic field, for different approximations of the above mentioned error covariance matrix. Also, 
in the context of a simplified one-dimensional MHD model, which retains part of the ingredients 
that make the complexity (and the beauty) of the geodynamo, Sun et al. (20071 have applied 
an optimal interpolation scheme that uses a Monte-Carlo method to calculate the same matrix, 
and studied the response of the system to assimilation for different temporal and spatial sampling 
frequencies. Both studies show a positive response of the system to SDA (i.e. better forecasts). 

In our opinion, though, SDA is strongly penalized by its formal impossibility to use current 
observations to improve past data records -even if this does not hamper its potential to produce 
good estimates of future core states. As said above, most of the information we have about the 
core is less that 500 yr old ( Jackson et al. 2000 1 . This record contains the signatures of the 



phenomena responsible for its short-term dynamics, possibly hydromagnetic waves with periods 
of several tens of years (Finlay and Jackson 2003 1. Our goal is to explore the VDA route in order 
to see to which extent high-resolution satellite measurements of the earth's magnetic field can 
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help improve the historical magnetic database, and identify more precisely physical phenomena 
responsible for short-term geomagnetic variations. To tackle this problem, we need a dynamical 
model of the high-frequency dynamics of the core, and an assimilation strategy. The aim of this 
paper is to reveal the latter, and illustrate it with a simplified one-dimensional nonlinear MHD 
model. Such a toy model, similar to the one used by Sun et al. (20071, retains part of the physics, 
at the benefit of a negligible computational cost. It enables intensive testing of the assimilation 
algorithm. 

This paper is organized as follows: the methodology we shall pursue in applying variational 
data assimilation to the geomagnetic secular variation is presented in Sect. [2] its implementation 
for the one-dimensional, nonlinear MHD toy model is described in detail in Sect. [3] Various 
synthetic assimilation experiments are presented in Sect . |4] the results of which are summarized 
and further discussed in Sect. H] 



2 Methodology 



In this section, we outline the bases of variational geomagnetic data assimilation, with the mid- 
term intent of improving the quality of the past geomagnetic record using the high-resolution 
information recorded by satellites. We resort to the unified set of notations proposed by Ide 
et al. (19971. What follows is essentially a transcription of the landmark paper by^Talagrand and 



Courtier (19871 with these conventions, transcription to which we add the possibility of imposing 



constraints to the core state itself during the assimilation process. 



2.1 Forward model 

Assume we have a prognostic, nonlinear, numerical model M which describes the dynamical 
evolution of the core state at any discrete time ti,i ^ {0, . . . , n}. If At denotes the time-step size, 
the width of the time window considered here is tn — to — nAt, the initial (final) time being to 
(tn). In formal assimilation parlance, this is written as 

Xi+i = Afi[xi], (1) 

in which x is a column vector describing the model state. If M relies for instance on the dis- 
cretization of the equations governing secular variation with a grid-based approach, this vector 
contains the values of all the field variables at every grid point. The secular variation equations 
could involve terms with a known, explicit time dependence, hence the dependence of M on time 
in Eq. ([T]). Within this framework, the modeled secular variation is entirely controlled by the 
initial state of the core, xq. 



2.2 Observations 



Assume now that we have knowledge of the true dynamical state of the core x* through databases 
of observations y° collected at discrete locations in space and time: 



(2) 



in which Hi and are the discrete observation operator and noise, respectively. For GSV, ob- 
servations consist of (scalar or vector) measurements of the magnetic field, possibly supplemented 
by decadal timeseries of the length of day, since these are related to the angular momentum of 
the core (Jault et al. 1988 Bloxham 1998). The observation operator is assumed linear and 
time-dependent: in the context of geomagnetic data assimilation, we can safely anticipate that its 
dimension will increase dramatically when entering the recent satellite era (1999-present). How- 
ever, H will always produce vectors whose dimension is much lower than the dimension of the 
state itself: this fundamental problem of undersampling is at the heart of the development of data 
assimilation strategies. The observational error is time-dependent as well: it is assumed to have 
zero mean and we denote its covariance matrix at discrete time ti by R^. 
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2.3 Quadratic misfit functions 

Variational assimilation aims here at improving the definition of the initial state of the core xq 
to produce modeled observations as close as possible to the observations of the true state. The 
distance between observations and predictions is measured using a quadratic misfit function Jh 

n 

Jh = Y. [^'^« - y°r Rr' [^.x, - y°] , (3) 

in which the superscript 'T' means transpose. In addition to the distance between observations 
and predictions of the past record, we might as well wish to try and apply some further constraints 
on the core state that we seek, through the addition of an extra cost function Jc 

n 

Jc = ^xfC7x„ (4) 

in which C is a matrix which describes the constraint one would like x to be subject to. This 
constraint can originate from some a priori ideas about the physics of the true state of the sys- 
tem, and its implication on the state itself, should this physics not be properly accounted for by 
the model M, most likely because of its computational cost. In the context of geomagnetic data 
assimilation, this a priori constraint can come for example from the assumption that fluid motions 
inside the rapidly rotating core are almost invariant along the direction of earth's rotation, ac- 



cording to Taylor-Proudman's theorem (e.g. Greenspan 1990). We shall provide the reader with 



an example for C when applying these theoretical concepts to the ID MHD model (see Sect. 4.2 1. 
Consequently, we write the total misfit function J as 

J^'^Jh+'^Jc. (5) 

where an and uq are the weights of the observational and constraint-based misfits, respectively. 
These two coefficients should be normalized; we will discuss the normalization in Sect. |4] 

2.4 Sensitivity to the initial conditions 

To minimize J, we express its sensitivity to xq, namely Vxq J. With our conventions, Vxo^/^ is a 
row vector, since a change in xq, (5xo, is responsible for a change in J, 5J , given by 

W=Vx„J-5xo. (6) 

To compute this gradient, we first introduce the tangent linear operator which relates a change 
in Xj_|_i to a change in the core state at the preceding discrete time, x^: 

(5x,+i = M[5Tii. (7) 

The tangent linear operator M[ is obtained by linearizing the model Mi about the state x^. 
Successive applications of the above relationship allow us to relate perturbations of the state 
vector Xi at a given model time ti to perturbations of the initial state xg: 

= Y[ A/j(5xo, e {1, . . . , n} (8) 

3=0 

The sensitivity of J to any x^ expresses itself via 

W=Vx,J-<5x„ (9) 

that is 

SJ=V^J-1[m'^5^o, {!,..., n}. (10) 

3=0 
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Additionally, after differentiating Eq. Q using Eqs. ^ and Q, we obtain 

Vx, J = aH{H,ic, - y°fR-^H,, + ac^fC, i e {0, . . . , n}. 

Gathering the observational and constraint contributions to J originating from every state vector 
Xi finally yields 

n i—1 

1 i— 1 



+aH(ffoXo - yafR-o ^Hq + acxj^C Uxq, 



which implies in turn that 



^ [aH(i/,x, - y'^fRi'H, + acxf C] J] Mj 

i=l j=0 

(i/oxo - yo)^Ro"^-ffo + acxj^C. 



(11) 



2.5 The adjoint model 



The computation of Vxo^ via Eq. (Ill is injected in an iterative method to adjust the initial state 



of the system to try and minimize J. The I + 1-th step of this algorithm is given in general terms 

by 

(12) 



X, 



i+i 



in which d is a descent direction, and an appropriate chosen scalar. In the case of the steepest 
descent algorithm, d' = (V^' J)^ , and is an a priori set constant. The descent direction is 
a column vector, hence the need to take the transpose of V^; J- In practice, the transpose of 



Eq. (|11|) yields, at the ^-th step of the algorithm, 

n 

= ■ ■ ■ ^^-1 [o^HHlRr\H,^\ - y°) + acC^^ 



V„i J 

■^0 



z=l 



(13) 



Introducing the adjoint variable a, the calculation of (V^i J)'^ is therefore performed practically 
by integrating the so-called adjoint model 



a^i = M'.-ia- + aHi/,^iR,_\(i/,_ixti - y°-i) + acCx^i, 
starting from ajj,,,]^ — 0, and going backwards in order to finally estimate 

(Vx<J)^-a[,. 



(14) 



(15) 



Equation (14 1 is at the heart of variational data assimilation (Talagrand 19971. Some remarks 



and comments concerning this so-called adjoint equation are in order: 

1. It requires to implement the transpose of the tangent linear operator, the so-called adjoint 
operator, M'J . If the discretized forward model is cast in terms of matrix-matrix and/or 
matrix- vector products, then this implementation can be rather straightforward (see Sect.js]). 
Still, for realistic applications, deriving the discrete adjoint equation can be rather convoluted 



(e.g. Bennett 2002 Chap. 4) 
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The discrete adjoint equation (Eq. 14 1 is based on the already discretized model of the 
secular variation. Such an approach is essentially motivated by practical reasons, assuming 
that we already have a numerical model of the geomagnetic secular variation at hand. We 
should mention here that a similar effort can be performed at the continuous level, before 
discretization. The misfit can be defined at this level; the calculus of its variations gives then 
rise to the Euler-Lagrange equations, one of which being the continuous backward, or adjoint, 
equation. One could then simply discretize this equation, using the same numerical approach 
as the one used for the forward model, and use this tool to adjust xq. According to [Bennett] 



(20021, though, the "discrete adjoint equation" is not the "adjoint discrete equation", the 



former breaking adjoint symmetry, which results in a solution being suboptimal (Bennett 



2002 § 4.1.6). 



3. Aside from the initial state Xq, one can in principle add model parameters (p, say) as 
adjustable variables, and invert jointly for xq and p, at the expense of expressing the discrete 
sensitivity of J to p as well. For geomagnetic VDA, this versatility might be of interest, in 
order for instance to assess the importance of magnetic diffusion over the time window of 
the historical geomagnetic record. 

4. The whole sequence of core states x' , i S {0, . . . ,n}, has to be kept in memory. This memory 
requirement can become quite significant when considering dynamical models of the GSV. 
Besides, even if the computational cost of the adjoint model is by construction equivalent 
to the cost of the forward model, the variational assimilation algorithm presented here is 
at least one or two orders of magnitude more expensive than a single forward realization, 
because of the number of iterations needed to obtain a significant reduction of the misfit 
function. When tackling 'real' problems in the future (as opposed to the illustrative problem 
of the next sections), memory and CPU time constraints might make it necessary to lower 
the resolution of the forward (and adjoint) models, by taking parameters values further 
away from the real core. A constraint such as the one imposed through Eq. ^ can then 
appear as a way to ease the pain and not to sacrifice too much physics, at negligible extra 
computational cost. 

We give a practical illustration of these ideas and concepts in the next two sections. 



3 Application to a one-dimensional nonlinear MHD model 

We consider a conducting fluid, whose state is fully characterized by two scalar flelds, u and b. 
Formally, b represents the magnetic field (it can be observed), and u is the velocity field (it is 
invisible) . 

3.1 The forward model 
3.1.1 Governing equations 

The conducting fluid has density p, kinematic viscosity v, electrical conductivity a, magnetic 
diffusivity 77, and magnetic permeability fi {rj = 1/ fia). Its pseudo-velocity u and pseudo-magnetic 
fleld b are both scalar flelds, deflned over a domain of length 2L, [—L, L]. We refer to pseudo flelds 
here since these flelds are not divergence-free. If they were so, they would have to be constant 
over the domain, which would considerably limit their interest from the assimilation standpoint. 
Bearing this remark in mind, we shall omit the 'pseudo' adjective in the remainder of this study. 

We choose L as the length scale, the magnetic diffusion time scale L'^/rj as the time scale, Bq 
as the magnetic fleld scale, and Bo/y/pfl as the velocity scale (i.e. the Alfven wave speed). Ac- 
cordingly, the evolution of u and b is controlled by the following set of non-dimensional equations: 

V(x,t) e]-l,l[x[0,T], 

dtu + S ud^u = S bd^b + Pmdlu, (16) 
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dtb + Sud^b = Sbd^u + dlb, 

supplemented by the boundary and initial conditions 

u{x, t) — if x = ±1, 
b{x,t) = ±lif2; = ±l, 

+ given = 0),6(-,< = 0). 



(17) 



(18) 
(19) 
(20) 



Eq. ( 16 1 is the momentum equation: the rate of change of the velocity is controlled by advection, 



magnetic forces and diffusion. Similarly, in the induction equation (17 1, the rate of change of the 
magnetic field results from the competition between inductive effects and ohmic diffusion. 
Two non-dimensional numbers define this system, 

which is the Lundquist number (ratio of the magnetic diffusion time scale to the Alfvcn time 
scale), and 

Pm ~ I'/ri, 

which is the magnetic Prandtl number, a material property very small for liquid metals - Pm ^ 10~^ 



for earth's core (e.g. Poirier 1988) 



3.1.2 Numerical model 

Fields are discretized in space using one Legendre spectral element of order N. In such a frame- 
work, basis functions are the Lagrangian interpolants /if defined over the collection of + 1 
Gauss-Lobatto-Legendre (GLL) points G {0,...,iV} (for a comprehensive description of 



the spectral element method, see Deville et al. 2002[ ). Figure [T] shows such a basis function for 
i = 50, N = 150. Having basis functions defined everywhere over [—1, 1] makes it straightforward 
to define numerically the observation operator H (see Sect. |3.3[ ). We now drop the superscript 
N for the sake of brevity. The semi-discretized velocity and magnetic fields are column vectors, 
denoted with bold fonts 



u(i) = [u{^o = -l,t),u{^i,t),...,u{^N = l,t)Y 
h{t) = [6(^0 = -l,t), 6(6, i),.-.,KCjv = l,i)]^ 



(21) 
(22) 



Discretization is performed in time with a semi-implicit finite-differencing scheme of order 1, 
explicit for nonlinear terms, and implicit for diffusive terms. As in the previous section, assuming 
that At is the time step size, we define ti = iAt,Ui — u(t — ti),hi = h{t = ti),i G {0, . . . ,n}. 
As a result of discretization in both space and time, the model is advanced in time by solving the 
following algebraic system 







■ 




■ f . " 






. H,-i _ 







where 



f 



= M/At + PmK, 

= M/A<+K, 

= M (ui/A< - 5ui Dui 

= M (bj/A< - Su, Db, 



5b, Db,) , 
5b, Du,) , 



(23) 



(24) 
(25) 
(26) 
(27) 



are the Helmholtz operators acting on velocity field and the magnetic field, and the forcing terms 
for each of these two, respectively. We have introduced the following definitions : 



which is the diagonal mass matrix. 
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Figure 1: An example of a basis function used to discretize the MHD model in space. This 
particular Lagrangian interpolant, hl^Q^, is obtained for a polynomial order N — 150, and it is 
attached to the 51st Gauss-Lobatto-Legendre point. 



• K, which is the so-called stiffness matrix (it is symmetric definite positive), 

• 0, which denotes the Hadamard product: (b u)^ = (u b)^ = bkUk, 



and D, the so-called derivative matrix 



dhl 



(28) 



the knowledge of which is required to evaluate the nonlinear terms. Advancing in time requires to 
invert both Helmholtz operators, which we do directly resorting to standard linear algebra routines 



(Anderson et al. 19991. Let us also bear in mind that the Helmholtz operators are symmetric (i.e. 
self-adjoint). 

In assimilation parlance, and according to the conventions introduced in the previous section, 
the state vector x is consequently equal to [u,b]-^, and its dimension is s = 2{N — 1) (since the 
value of both the velocity and magnetic fields are prescribed on the boundaries of the domain) . 



3.2 The true state 

Since we are dealing in this paper with synthetic observations, it is necessary to define the true 
state of the ID system as the state obtained via the integration of the numerical model defined in 
the preceding paragraph, for a given set of initial conditions, and specific values of the Lundquist 
and magnetic Prandtl numbers, S and Pm. The true state (denoted with the superscript 't') will 
always refer to the following initial conditions 

u*{x, t^O) ^ sm{nx) + (2/5) sin(57ra;), (29) 
b\x,t = 0) = cos(7ra;) -H2sin[7r(a;-|- 1)/4], (30) 

along with S — 1 and Pm — 10^''. The model is integrated forward in time until T — 0.2 (a fifth 
of a magnetic diffusion time). The polynomial order used to compute the true state is iV = 300, 
and the time step size At = 2 10~^. Figure [2] shows the velocity (left) and magnetic field (right) 
at initial (black curves) and final (red curves) model times. The low value of the magnetic Prandtl 
number Pm reflects itself in the sharp velocity boundary layers that develop near the domain 
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Figure 2: The true state used for synthetic variational assimilation experiments. Left: the first, 
t = (black) and last, t = T (red) velocity fields. Right: the first, t = (black) and last, t = T 
(red) magnetic fields. 



boundaries, while the magnetic field exhibits in contrast a smooth profile (the magnetic diffusivity 
being three orders of magnitude larger than the kinematic viscosity). To properly resolve these 
Hartmann boundary layers there must be enough points in the vicinity of the domain boundaries: 
we benefit here from the clustering of GLL points near the boundaries ( Deville et al.[ 20021 



Besides, even if the magnetic profile is very smooth, one can nevertheless point out here and there 
kinks in the final profile. These kinks are associated with sharp velocity gradients (such as the 



one around x = 0.75) and are a consequence of the nonlinear bd^u term in the induction Eq. (17 1 



3.3 Observation of the true state 

In order to mimic the situation relevant for the earth's core and geomagnetic secular variation, 
assume that we have knowledge of b at discrete locations in space and time, and that the velocity 
u is not measurable. For the sake of generality, observations of b are not necessarily made at 
collocation points, hence the need to define a spatial observation operator Hf (at discrete time 
U) consistent with the numerical approximation introduced above. If nf denotes the number of 
virtual magnetic stations at time ti, and their locations {j £{!,..., nf}), H~ is a rectangular 
X {N + 1) matrix, whose coefficients write 

Hli = h^{£.l,). (31) 

A database of magnetic observations y° = b° is therefore produced at discrete time ti via the 
matrix-vector product 

b° = Hfh\. (32) 
Integration of the adjoint model also requires the knowledge of the transpose of the observation 



operator (Eq. 14 1 , the construction of which is straightforward according to the previous definition. 
To construct the set of synthetic observations, we take for simplicity the observational noise to 
be zero. During the assimilation process, we shall assume that estimation errors are uncorrelated, 
and that the level of confidence is the same for each virtual observatory. Consequently, 

R» = 1°, (33) 

in which 1° is the x identity matrix, throughout the numerical experiments. 

As an aside, let us notice that magnetic observations could equivalently consist of an (arbitrarily 
truncated) set of spectral coefficients, resulting from the expansion of the magnetic field on the 
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basis of Legendre polynomials. Our use of stations is essentially motivated by the fact that our 
forward model is built in physical space. For real applications, a spectral approach is interesting 
since it can naturally account for the potential character of the field in a source-free region; 
however, it is less amenable to the spatial description of observation errors, if these do not vary 
smoothly. 



3.4 The adjoint model 

3.4.1 The tangent linear operator 

As stated in the the previous section, the tangent linear operator M^' at discrete time ti is obtained 
at the discrete level by linearizing the model about the current solution (u^jb.;). By perturbing 
these two fields 



b, 



^b, 



5u,, 
- Shi, 



(34) 
(35) 



we get (after some algebra) 



having introduced the (A^ + 1)^ following matrices 



<5ui+i 





















A, 

E,; 



(l/A<-S'Du, ©-S-u, 0D), 
{Sh, D + 5Db,0) , 
{~SDh, -5b, D) , 
(I/At - Su, D + SDu,Q) . 



(36) 
(37) 
(38) 
(39) 



Aside from the (A'^ + 1)^ identity matrix I, matrices and notations appearing in these definitions 
have already been introduced in |3.1.2[ In connection with the general definition introduced in 
the previous section, 5x^+1 = Al^Sxi, M- is the block matrix 



A. 
Q 



B, 
E, 



(40) 



3.4.2 Implementation of the adjoint equation 



The sensitivity of the model to its initial conditions is computed by applying the adjoint operator, 
Mf", to the adjoint variables - see Eq. (14 1. According to Eq. (40), one gets 



M: 



IT 



Af 
BF 



E^ 



with each transpose given by 

Af 
Bf 



(l/Ai - 5u, D"^ - 5D^u,0) MH-i, 
(5D^b, +5b, D^) MH;1, 
(-5b, - 5D^b,0) MH-\ 
{\/M - 5D'^Ui +5u, D^) MH(7\ 



(41) 



(42) 
(43) 
(44) 
(45) 



In writing the equation in this form, we have used the symmetry properties of the Helmholtz 
and mass matrices, and introduced the transpose of the derivative matrix, Z)-^. Programming the 
adjoint model is very similar to programming the forward model, provided that one has cast the 
latter in terms of matrix-matrix, matrix- vector, and Hadamard products. 
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4 Synthetic assimilation experiments 



Having all the numerical tools at hand, we start out by assuming that we have imperfect knowledge 
of the initial model state, through an initial guess Xg, with the model parameters and resolution 
equal to the ones that helped us define the true state of |3.2| We wish here to quantify how 
assimilation of observations can help improve the knowledge of the initial (and subsequent) states, 
with particular emphasis on the influence of spatial and temporal sampling. In the series of results 
reported in this paragraph, the initial guess at model initial time is : 



u^{x,t = 0) = sin(7ra;), 

b\x,t^O) = cos(7r2:) +2sin[7r(2: + l)/4] + (l/2)sin(27rx). 



(46) 
(47) 



With respect to the true state at the initial time, the first guess is missing the small-scale compo- 
nent of u, i.e. the second term on the right-hand side of Eq. (29 1. In addition, our estimate of 6 



has an extra parasitic large-scale component (the third term on the right-hand side of Eq. (47 1), 
a situation that could occur when dealing with the GSV, for which the importance of unmodeled 



small-scale features has been recently put forward given the accuracy of satellite data (Eymin 
and Hulot , 2005 1 . Figure [s] shows the initial and final and , along with u* and 6* at the 



same epochs for comparison, and the difference between the two, multiplied by a factor of five. 
Differences in b are not pronounced. Over the time window considered here, the parasitic small- 
scale component has undergone considerable diffusion. To quantify the differences between the 
true state and the guess, we resort to the L2 norm 



11/11 = 




Pdx, 



and define the relative magnetic and fluid errors at time ti by 



t 



/ 11^*11. 



(48) 
(49) 



The initial guess given by Eqs. (46 1(47 1 is characterized by the following errors: Cg — 21.6%, — 
2.9%, e^f = 37.1%, and el = 37.1%. 



4.1 Improvement of the initial guess with no a priori constraint on the 
state 

4.1.1 Regular space and time sampling 

Observations of are performed at virtual observatories which are equidistant in space, at a 
number of epochs nt evenly distributed over the time interval. Assuming no a priori constraint on 
the state, we set ac = in Def. [sj The other constant an = l/{ntv^). 

The minimization problem is tackled by means of a conjugate gradient algorithm, a la Polak- 



Ribiere (Shewchuk 19941. Iterations are stopped either when the initial misflt has decreased by 8 
orders of magnitude, or when the iteration count exceeds 5,000. In most cases, the latter situation 
has appeared in our simulations. A typical minimization is characterized by a fast decrease in 
the misflt during the flrst few tens of iterations, followed by a slowly decreasing (almost flat) 
behaviour. Even if the solution keeps on getting better (i.e. closer to the synthetic reality) 
during this slow convergence period, practical considerations (having in mind the cost of future 
geomagnetic assimilations) prompted us to stop the minimization. 

A typical example of a variational assimilation result is shown in Fig. [4] In this case, = 20 
and nt = 20. The recovery of the flnal magnetic fleld 6„ is excellent (see Fig. l4ji), the relative 



L2 error being 1.8 10 . The beneflt here is double: flrst, the network of observatories is dense 
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Figure 3: Initial guesses used for the variational assimilation experiments, plotted against the 
corresponding true state variables. Also plotted is five times the difference between the two. a: 
velocity at time 0. b: velocity at final time T. c: magnetic field at time 0. d: magnetic field at 
final time T. In each panel, the true state is plotted with the black line, the guess with the green 
line, and the magnified difference with the blue line. 
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Figure 4: Synthetic assimilation results, a): velocity at initial model time t — 0. b): velocity at 
final time t = T. c): magnetic field at initial time t — 0. d): magnetic field at final time t = T. In 
each panel, the true field is plotted in black, the assimilated field (starting from the guess shown 
in Fig. ^ in green, and the difference between the two, multiplied by a factor of 5, is shown in 
blue. The red triangles indicate the location of the n'^ virtual magnetic observatories (n'^ = 20 in 
this particular case). 
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Figure 5: Dynamical evolution of L2 errors (logarithmic value) for the magnetic field (a) and the 
fluid velocity (b). Dashed lines : errors for initial guesses. Solid lines : errors after variational as- 
similation. Circles represent instants are which magnetic observations are made. In this particular 



case, nt — 20 and 



20. 



enough to sample properly the field, and second, a measurement is made exactly at this discrete 
time instant, leaving no time for error fields to develop. When the latter is possible, the recovered 
fields can be contaminated by small-scale features, that is features that have length scales smaller 
than the spatial sampling scale. We see this happening in Fig. [4j;), in which the magnified 
difference between the recovered and true bo, shown in blue, appears indeed quite spiky; Cq has 
still decreased from an initial value of 21.6% (Fig. |3j;) down to 1.2%. Results for velocity are 
shown in Figs. and|4j3. The recovered velocity is closer to the true state than the initial guess: 
this is the expected benefit from the nonlinear coupling between velocity and magnetic field in 
Eqs. (16)-pT|. The indirect knowledge we have of u, through the observation of b, is sufficient to 
get better estimates of this field variable. At the end of the assimilation process, Eq and e", which 
were approximately equal to 37% with the initial guess, have been brought down to 8.2 and 4.7 
%, respectively. The velocity at present time (Fig. |4]) is remarkably close to the true velocity, save 
for the left boundary layer sharp structure, which is undersampled (see the distribution of red 
triangles). We further document the dynamical evolution of L2 errors by plotting on Fig. [s] the 
temporal evolution of e'' and e" for this particular configuration. Instants at which observations 
are made are represented by circles, and the temporal evolution of the guess errors are also shown 
for comparison. The guess for the initial magnetic field is characterized by a decrease of the error 
that occurs over « .1 diffusion time scale, that is roughly the time it takes for most of the parasitic 
small-scale error component to diffuse away, the error being then dominated at later epochs by 
advection errors, originating from errors in the velocity field. The recovered magnetic field (Fig. [5^, 
solid line), is in very good agreement with the true field as soon as measurements are available 
{t > 1% of a magnetic diffusion time, see the circles on Fig. [5^). Even though no measurements 
are available for the initial epoch, the initial field has also been corrected significantly, as discussed 
above. In the latter parts of the record, oscillations in the magnetic error field are present -they 
disappear if the minimization is pushed further (not shown). 

The unobserved velocity field does not exhibit such a drastic reduction in error as soon as 
observations are available (Fig. [sja, solid line). Still, it is worth noticing that the velocity error is 
significantly smaller in the second part of the record, in connection with the physical observation 
that most of the parasitic small-scale component of the field has decayed away (see above): advec- 
tion errors dominate in determining the time derivative of b in Eq. ( 17l, leaving room for a better 
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Figure 6: Systematic study of the response of the one-dimensional MHD system to variational 
assimilation. Shown are the logarithms of L2 errors for the t — Q {a) and t — T {b) magnetic field, 
and the t — Q {c) and t — T {d) velocity field, versus the number of times observations are made 
over [0,T], ni, using spatial networks of varying density (n^ = 5, 10,20,50 and 100). 



assessment of the value of u. For other cases (different rP and nt), we find a similar behaviour 
(not shown). We comment on the effects of an irregular time sampling on the above observations 
in section 14.1.31 

— '. . . . . Q . 

Having in mind what one gets in this particular {nt,n ) configuration, we now summarize 

in Fig. [6] results obtained by varying systematically these 2 parameters. After assimilation, the 
logarithmic value of the L2 velocity and magnetic field errors, at the initial and final stages {i — 
and i — n), are plotted versus nt, using = 5, 10,20,50, and 100 virtual magnetic stations. As 
far as temporal sampling is concerned, nt can be equal to 1 (having one observation at present 
time only), 10, 20, 50 or 100. Inspection of Fig. [6] leads us to make the following comments: 

• Regarding h : 

— 50 stations are enough to properly sample the magnetic field in space. In this case 
Tit = 1 is sufficient to properly determine b„, and no real improvement is made when 
increasing nt (Fig. |6]d). During the iterative process, the value of the field is brought 
to its observed value at every station of the dense network, and this is it: no dynamical 
information is needed. 

— When, on the other hand, spatial sampling is not good enough, information on the 
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dynamics of b helps improve partially its knowledge at present time. For instance, we 
get a factor of 5 reduction in with = 20, going from nt — 1 to nt = 10 (Fig. 
circles). The improvement then stabilizes upon increasing nt: spatial error dominates. 

— This also applies for the initial magnetic field bp, see Fig. [6^. As a matter of fact, 
having no dynamical information about b (nt = 1) precludes any improvement on bo, 
for any density of the spatial network. Improvement occurs for nt > 1. If the spatial 
coverage is good enough (n^ > 50), no plateau is reached, since the agreement between 
the assimilated and true fields keeps on getting better, as it should. 

• Regarding u : 

— The recovered u is always sensitive to spatial resolution, even for nt — 1 (Figs, [ot and 

— If nt is increased, the error decreases and reaches a plateau which is again determined 
by spatial resolution. This holds for Bq and eJJ. For the reason stated above, u„ is 
better known than Uq. The error is dominated in both cases by a poor description of 
the left boundary layer (see the blue curves in Figs. [4^ and|4]3). The gradient associated 
with this layer is not sufficiently well constrained by magnetic observations (one reason 
being that the magnetic diffusivity is three times larger than the kinematic viscosity). 
Consequently, we can speculate that the error made in this specific region at the final 
time is retro-propagated and amplified going backwards in time, through the adjoint 
equation, resulting in Cq > e^. 

4.1.2 Irregular spatial sampling 

We have also studied the effect of an irregular spatial sampling by performing a suite of simulations 
identical to the ones described above, save that we assumed that stations were only located in the 
left half of the domain (i.e. the [—1,0] segment). 

The global conclusion is then the following: assimilation results in an improvement of estimates 
of b and u in the sampled region, whereas no benefit is visible in the unconstrained region. To 
illustrate this tendency (and keep a long story short), we only report in Fig. [7]the recovered u 
and b for {n^,nt) = (10, 20), which corresponds to the "regular" case depicted in Fig. [4] deprived 
from its 10 stations located in [0, 1]. The lack of transmission of information from the left-hand 
side of the domain to its right-hand side is related to the short duration of model integration 
(0.2 magnetic diffusion time, which corresponds to 0.2 advective diffusion time with our choice 
of 5 = 1). We shall comment further on the relevance of this remark for the assimilation of the 
historical geomagnetic secular variation in the discussion. 

The lack of observational constraint on the right-hand side of the domain results sometimes in 
final errors larger than the initial ones (compare in particular Figs. |7k andlTb, with Figs. Ek and 

We also note large error oscillations located at the interface between the left (sampled) and 
right (not sampled) regions, particularly at initial model time (Figs.|7^ and[7|:). The contrast in 
spatial coverage is likely to be the cause of these oscillations (for which we do not have a formal 
explanation); this type of behaviour should be kept in mind for future geomagnetic applications. 

4.1.3 Irregular time sampling 

We can also assume that the temporal sampling rate is not constant (keeping the spatial network of 
observatories homogeneous), restricting for instance drastically the epochs at which observations 
are made to the last 10 % of model integration time, the sampling rate being ideal (that is 
performing observations at each model step). Not surprisingly, we are penalized by our total 
ignorance of the 90 remaining per cent of the record. We illustrate the results obtained after 
assimilation with our now well-known array of n = 20 stations by plotting the evolution of the 
errors in b and u (as defined above) versus time in Fig.ls] Although the same amount of information 
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Figure 7: Synthetic assimilation results obtained with an asymmetric network of virtual obser- 
vatories (red triangles). Other model and assimilation parameters as in Fig. |4] a): velocity at 
initial model time t = 0. b): velocity at final model time t — T. c): magnetic field at t = 0. d): 
magnetic field at i = T. In each panel, the true field is plotted in black, the assimilated field in 
green, and the difference between the two, multiplied by a factor of 5, is shown in blue. 
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Figure 8: Same as Fig. [5] save that the nt = 20 epochs at which measurements are made are 
concentrated over the last 10% of model integration time. 

(n^nt — 400) has been collected to produce Fig s. H and [s] the uneven temporal sampling of the 
latter has dramatic consequences on the improvement of the estimate of h. In particular, the initial 
error eg remains large. The error decreases then linearly with time until the first measurement is 
made. We also observe that the minimum is obtained in the middle of the observation era. The 
poor quality of the temporal sampling, coupled with the not-sufficient spatial resolution obtained 
with these 20 stations, does not allow us to reach error levels as small as the ones obtained in Fig. 
[5] even at epochs during which observations are made. The velocity is sensitive to a lesser extent 
to this effect, with velocity errors being roughly 2 times larger in Fig. [8] than in Fig. [5] 



4.2 Imposing an a priori constraint on the state 

As stated in Sect.|2] future applications of variational data assimilation to the geomagnetic secular 
variation might require to try and impose a priori constraints on the core state. In a kinematic 
framework, this is currently done in order to restrict the extent of the null space when trying to 



invert for the core flow responsible for the GSV (Backus 1968 Le Mouel 19841 



Assume for instance that we want to try and minimize the gradients of the velocity and 
magnetic fields, in a proportion given by the ratio of their diffusivities, that is the magnetic 
Prandtl number Pm, at any model time. The associated cost function is written 



i=0 



[bf D^Db, + Pm (uf D^Dui)] 



(50) 



in which D is the derivative matrix introduced in §3.1.2[ The total misfit reads, according to 
Eq. (|5| 

J = anJu + acJc, 

with aff — l/{ntn ) as before, and ac = (3/[n{N — 1)], in which (3 is the parameter that controls 
the constraint to observation weight ratio. Response of the assimilated model to the imposed 
constraint is illustrated in Fig. [9] using the {nt — 20, = 20) reference case of Fig. [4] for three 
increasing values of the P parameter: 10~^, 1, and 10^, and showing also for reference what happens 
when /3 = 0. We show the error fields (the scale is arbitrary, but the same for all curves) at the 
initial model time, for velocity (left panel) and magnetic field (right panel). The L2 errors for each 
field at the end of assimilation indicate that this particular constraint can result in marginally 
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Figure 9: Influence of an a priori imposed constraint (in this case aiming at reducing the gradients 
in the model state) on the results of variational assimilation. Shown are the difference fields 
(arbitrary scales) between the assimilated and true states, for the velocity field (left panel) and 
the magnetic field (right panel), at initial model time. Again, as in Fig. Ill we have made nt — 
20 measurements at n'^ — 20 evenly distributed stations. /3 measures the relative ratio of the 
constraint to the observations. Indicated for reference are the L2 errors corresponding to each 
configuration. The grey line is the zero line. 
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iteration count 

Figure 10: Convergence behaviour for different constraint levels (3. The ratio of the current value 
of the misfit j' (normalized by its initial value J°) is plotted against the iteration count /. (3 
measures the strength of the constraint imposed on the state relative to the observations. 



better estimate of the initial state of the model, provided that the value of the parameter (3 is kept 
small. For (3 = 10""'^, the initial magnetic field is much smoother than the one obtained without 
the constraint and makes more physical sense (Fig. [9|i). The associated velocity field remains 
spiky, with peak to peak error amplitudes strongly reduced in the heart of the computational 
domain (Fig. [9]:). This results in smaller errors (reduction of about 20% for 6o and 10% for uq). 
Increasing further the value of (3 leads to a magnetic field that is too smooth (and an error field 
even dominated by large-scale oscillations, see Fig. [9]i), simply because too much weight has been 
put on the large-scale components of b. The velocity error is now also smooth (Fig. [9^), at the 
expense of a velocity field being further away from the sought solution (eg = 11.7%), especially 
in the left Hartmann boundary layer. In the case of real data assimilation (as opposed to the 
synthetic case here, the true state of which we know, and departures from which we can easily 
quantify), we do not know the true state. To get a feeling for the response of the system to the 
imposition of an extra constraint, it is nevertheless possible to monitor for instance the convergence 
behaviour during the descent. On Fig.jlOj the ratio of the misfit to its initial value is plotted versus 
the iteration number in the conjugate gradient algorithm (log-log plot). If [3 is small, the misfit 
keeps on decreasing, even after 5,000 iterations (green curve). On the other hand, a too strong a 
constraint (blue and red curves in Fig. 10 1 is not well accommodated by the model and results in 
a rapid fiattening of the convergence curve, showing that convergence behaviour can be used as a 
proxy to assess the efficacy of an a priori imposed constraint. 

Again, we have used the constraint given by Eq. (50 1 for illustrative purposes, and do not claim 
that this specific low-pass filter is mandatory for the assimilation of GSV data. Similar types of 
constraints are used to solve the kinematic inverse problem of GSV (Bloxham and Jackson 1991 1; 



see also Pais et al. (2004) and Amit and Olson (2004| for recent innovative studies on the subject. 
The example developed in this section aims at showing that a formal continuity exists between 
the kinematic and dynamical approaches to the GSV. 
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4.3 Convergence issues 



In most of the cases presented above, the iteration counts had reached 5,000 before the cost 
function had decreased by 8 orders of magnitude. Even though the aim of this paper is not to 
address specifically the matter of convergence acceleration algorithms, a few comments are in 
order, since 5,000 is too large a number when considering two- or three-dimensional applications. 

• In many cases, a reduction of the initial misfit by only 4 orders of magnitude gives rise to 
decent solutions, obtained typically in a few hundreds of iterations. For example, in the case 
corresponding to Fig. |4] a decrease of the initial misfit by 4 orders of magnitude is obtained 
after 475 iterations. The resulting error levels are already acceptable : — 12 10^^, 

= 7.5 10^2^ e|5 = 1.8 lO^^ and = 3.0 lO"''. 

• More importantly, in future applications, convergence will be sped up through the intro- 
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duction of a background error covariance matrix B, resulting in an extra term (Ide et al. 

^[xo - Xb]^B"^[xo - Xf,] 



added to the cost function (Eq. ([5|). Here, Xf, denotes the background state at model time 0, 
the definition of which depends on the problem of interest. In order to illustrate how this 
extra term can accelerate the inversion process, we have performed the following assimilation 
experiment: we take the network of virtual observatories of Fig. [4] and define the background 
state at model time to be zero for the velocity field (which is not directly measured) , and the 
polynomial extrapolation of the t — magnetic observations made at the n — 20 stations 
on the N + 1 GLL grid points for the magnetic field (resorting to Lagrangian interpolants 
defined by the network of stations). The background error covariance matrix is chosen to 
be diagonal, without cross-covariance terms. This approach enables a misfit reduction by 
5 orders of magnitude in 238 iterations, with the following L2 error levels : Cq = 13 10^^, 
= 11.9 10~^, eg = 2.6 10^^, and = 2.6 10^"*. This rather crude approach is beneficial 
for a) the computational cost and b) the estimate of the magnetic field. The recovery of 
the velocity is not as good as it should be, because we have made no assumption at all on 
the background velocity field. In future applications of VDA to the GSV, some a priori 
information on the background velocity field inside the core will have to be introduced in 
the assimilation process. The exact nature of this information is beyond the scope of this 
study. 



5 Summary and conclusion 

We have laid the theoretical and technical bases necessary to apply variational data assimilation 
to the geomagnetic secular variation, with the intent of improving the quality of the historical geo- 
magnetic record. For the purpose of illustration, we have adapted these concepts (well established 
in the oceanographic and atmospheric communities) to a one-dimensional nonlinear MHD model. 
Leaving aside the technical details exposed in section [3] we can summarize our findings and ideas 
for future developments as follows: 

• Observations of the magnetic field always have a positive impact on the estimate of the 
invisible velocity field, even if these two fields live at different length scales (as could be 
expected from the small value of the magnetic Prandtl number). 

• With respect to a purely kinematic approach, having successive observations dynamically 
related by the model allows one to partially overcome errors due to a poor spatial sampling 
of the magnetic field. This is particularly encouraging in the prospect of assimilating main 
geomagnetic field data, the resolution of which is limited to spherical harmonic degree 14 
(say), because of (remanent or induced) crustal magnetization. 
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Figure 11: Dynamical evolution of L2 errors (logarithmic value) for the magnetic field (a) and the 
fluid velocity (b). Black lines : errors for initial guesses. Green (red) lines : errors for assimilation 
results that do (not) incorporate the data obtained by a dense virtual network of magnetic stations, 
which aims at mimicking the satellite era -the blue segment on each panel-, spanning the last 5 % 
of model integration time. 



Over the model integration time (20 % of an advection time), regions poorly covered ex- 
hibit poor recoveries of the true fields, since information does not have enough time to be 
transported there from well covered regions. In this respect, model dynamics clearly controls 
assimilation behaviour. Concerning the true GSV, the time window we referred to in the 
introduction has a width of roughly a quarter of an advective time scale. Again, this is rather 
short to circumvent the spatial limitations mentioned above, if advective transport controls 
the GSV catalog. This catalog, however, could contain the signature of global hydromag- 



netic oscillations (Hide 1966 Finlay and Jackson 20031, in which case our hope is that 



problems due to short duration and coarse spatial sampling should be alleviated. This issue 
is currently under investigation in our simplified framework, since the toy model presented 
here supports Alfven waves. 



• A priori imposed constraints (such as the low-pass filter of Sect. 4.2 1 can improve assimilation 
results. They make variational data assimilation appear in the formal continuity of kinematic 
geomagnetic inverse problems as addressed by the geomagnetic community over the past 40 
years. 

Finally, in order to illustrate the potential interest of applying VDA techniques to try and improve 



the recent GSV record, we show in Fig. 11 the results of two synthetic assimilation experiments. 
These are analogous to the ones described in great length in Sect. |4] (same physical and numerical 
parameters, constraint parameter (3 = 10~^). In both cases, observations are made by a network 
of 6 evenly distributed stations during the first half of model integration time (the logbooks era, 
say) . The second half of the record is then produced by a network of 15 stations for case A (the 
observatory era). For case B, this is also the case, save that the last 5% of the record are obtained 
via a high-resolution network of 60 stations. The two records therefore only differ in the last 5% 
of model integration time. Case B is meant to estimate the potential impact of the recent satellite 
era on our description of the historical record. 

The evolution of the magnetic error backwards in time (Fig. Ill) shows that the benefit 
due to the dense network is noticeable over three quarters of model integration time, with an 
error reduction of roughly a factor of 5. The velocity field is (as usual) less sensitive to the better 
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quality of the record; still, it responds well to it, with an average decrease of e" on the order of 
20%, evenly distributed over the time window. 

Even if obtained with a simplified model (bearing in particular in mind that real geomagnetic 
observations are only available at the core surface), these results are promising and indicate that 
VDA should certainly be considered as the natural way of using high-quality satellite data to 



refine the historical geomagnetic record in order to 'reassimilate' (Talagrand 19971 pre-satellite 
observations. To do so, a good initial guess is needed, which is already available (Jackson et al. 



2000); also required is a forward model (and its adjoint) describing the high-frequency physics 
of the core. This model could either be a full three-dimensional model of the geodynamo, or a 
two-dimensional, specific model of short-period core dynamics, based on the assumption that this 
dynamics is quasi-geostrophic (Jault 20061. The latter possibility is under investigation. 
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